/*
 * aln.c
 *
 *  Created on: Oct 25, 2013
 *      Author: jiang
 */

#include <zlib.h>
#include <stdlib.h>
#include <stdio.h>
#include <time.h>

#include "kseq.h"
#include "util.h"
#include "bntseq.h"
#include "seqio.h"
#include "aln.h"
#include "utils.h"
#include "lookup.h"
#include "global.h"
#include "khash.h"

#include "stdaln.h"
#include "aln.h"

//DEBUG ONLY//////
#include <stdbool.h>

#include <limits.h>
/////////////////

#include <inttypes.h>

///DEBUG ONLY ////////
char *int2bin(unsigned n, char *buf) {
#define BITS (sizeof(n) * CHAR_BIT)

	static char static_buf[BITS + 1];
	int i;

	if (buf == NULL)
		buf = static_buf;

	for (i = BITS - 1; i >= 0; --i) {
		buf[i] = (n & 1) ? '1' : '0';
		n >>= 1;
	}

	buf[BITS] = '\0';
	return buf;

#undef BITS
} ////

KHASH_MAP_INIT_INT64(64, char)
//DEBUG ONLY///////////////////////////
uint32_t gloablCount = 0;
uint64_t NinReadsCount = 0;
uint32_t gloabl3KmerCount = 0;
uint32_t unableMapReadsCount = 0;
uint64_t debug_totalKmer = 0;
uint64_t scoreLowThresCount = 0;
uint64_t notHasDKmer = 0;

#define __gen_ap(par, opt) do {									\
		int i;													\
		for (i = 0; i < 25; ++i) (par).matrix[i] = -(opt)->b;	\
		for (i = 0; i < 4; ++i) (par).matrix[i*5+i] = (opt)->a; \
		(par).gap_open = (opt)->q; (par).gap_ext = (opt)->r;	\
		(par).gap_end = (opt)->r;								\
		(par).row = 5; (par).band_width = opt->bw;				\
	} while (0)

unsigned char d2char_table[4] = { 'A', 'C', 'G', 'T' };

static void bwaReadsWrite(FILE *bwaReadsFile, read_t *read) {

	int i;
	fseek(bwaReadsFile, 0, SEEK_END);
	fprintf(bwaReadsFile, "@%s \n", read->name);
	for (i = 0; i < read->len; i++) {
		fprintf(bwaReadsFile, "%c", d2char_table[read->seq[i]]);
	}
	fprintf(bwaReadsFile, "\n+\n");
	for (i = 0; i < read->qual_len; i++) {
		fprintf(bwaReadsFile, "%c", read->qual[i]);
	}
	fprintf(bwaReadsFile, "\n");

}

static uint32_t add2UnmapedReads(FILE *bwaReadsFile, read_t *read) {
	unableMapReadsCount++;
	//For  these read, we use hash method.
	//First we get statistics that which region in the reference a
	bwaReadsWrite(bwaReadsFile, read);
	return 0;
}

static uint32_t lowScoreHandle(FILE *bwaReadsFile, read_t *read) {
	scoreLowThresCount++;
	bwaReadsWrite(bwaReadsFile, read);
	return 0;
}

static uint32_t tripN(ubyte_t ch) {
	srand(time(NULL));
	if (ch > 3) {
		NinReadsCount += 1;
		return rand() % 4;
	} else
		return ch;
}

static int64_t dmap_extend(uint32_t *len, const uint8_t *seq1,
		const uint8_t *seq2, const AlnParam *ap, int thres, int len1, int len2,
		opt_t *opt) {
	int matrix[25];
	AlnParam par;
	int subo, path_len, score, n_cigar, i;
	uint32_t * cigar;
	path_t *path;

	par.matrix = matrix;
	__gen_ap(par, opt);

	path = (path_t *) malloc((len1 + len2) * sizeof(path_t));

	score = aln_local_core(seq1, len1, seq2, len2, &par, path, &path_len, 1,
			&subo);
	cigar = aln_path2cigar32(path, path_len, &n_cigar);

	/*
	 fprintf(stderr,"LEN %d seq1:",len1);
	 for( i = 0; i < len1; i++)
	 fprintf(stderr,"%d",seq1[i]);
	 fprintf(stderr,"\n");
	 fprintf(stderr,"LEN %d seq2:",len2);
	 for( i = 0; i < len2; i++)
	 fprintf(stderr,"%d",seq2[i]);

	 fprintf(stderr,"\n");
	 for( i = 0; i < path_len; i++)
	 fprintf(stderr,"i,j (%d,%d) \n",path[i].i,path[i].j);

	 for( i = 0; i < n_cigar; i++)
	 fprintf(stderr,"%d  ", cigar[i]);

	 //fprintf(stderr,"score: %d  path_len: %d \n", score, path_len);

	 // subo means what??
	 /*
	 if (ret < 0 || subo == ret) { // no hit or tandem hits
	 free(path);
	 free(cigar);
	 free(seq1);
	 n_cigar = 0;
	 return 0;
	 }
	 */
	*len = path_len;

	free(path);
	return score;
}


static void printPos( uint64_t rawPos, int score,uint32_t readLength , int64_t l_ref, bntseq_t *bns,float mapQuality,uint64_t freq)
{
	int chrID, posSign;
	uint32_t ScorePos;
	posSign = (rawPos > BIGNUMBER) ? -1 : 1;

	ScorePos = (rawPos > BIGNUMBER) ? rawPos - BIGNUMBER :rawPos;

	if(posSign == -1){
		ScorePos = l_ref - ScorePos - readLength;
	}
	bns_cnt_ambi(bns, ScorePos, readLength, &chrID);
	ScorePos = ScorePos - bns->anns[chrID].offset + 1;
	//fprintf(stdout, "\t%d\t%llu\t%f\t%s", posSign, ScorePos,mapQuality,bns->anns[chrID].name);
	fprintf(stdout, " %llu_%d_%llu", ScorePos, score,freq);
}

/*
 * Do  alignment for a single read.
 */
static uint32_t Aln(bntseq_t *bns, alnParam parameter,
		pos_shift *positions_shifts, uint32_t dkmerN, const int64_t l_ref,
		const uint8_t *pac, const uint32_t readLength, const ubyte_t *read,
		opt_t *opt, misReadsFiles otherReadsFiles, read_t *readWhole) {

	int ret, i, distCount = 0, topZ, posFreqDicCount = 0;
	khiter_t k;
	uint8_t *left_target, *right_target,*mid_target;
	uint64_t distinctPositions[ARRAY_LEN];
	pair posFreqDic[ARRAY_LEN];
	posMatchedRegion_t posMatchedRegion[ARRAY_LEN]; // to present already matched region while seeding process
	uint32_t posMatchedRegionCount, mpos;

	int maxScore = -1000, maxScorePosSign = 0;
	uint32_t maxScorePos;
	float mapQuality, c1=0.5;
	int chrID,old_chaining = 0;
	int scores[ARRAY_LEN];


	//Get distinct positions
	khash_t(64) *h = kh_init(64);
	for (i = 0; i < dkmerN; i++) {
		k = kh_put(64, h, positions_shifts[i].position, &ret); // dictionary process
		kh_val(h,k) = 0;
	}
	for (i = 0; i < dkmerN; i++) {
		kh_val(h, kh_get(64, h, positions_shifts[i].position)) += 1;
	}

	//Pick top_Z high frequent
	for (k = kh_begin(h); k != kh_end(h); ++k) {
		if (kh_exist(h, k)) {
			distinctPositions[distCount++] = kh_key(h, k);
			posFreqDic[posFreqDicCount].first = kh_key(h, k);
			posFreqDic[posFreqDicCount].second = kh_val(h, k);
			posFreqDicCount++;
		}
	}

	kh_destroy(64, h);

	//Sort by frequency, descending order
	qsort(posFreqDic, posFreqDicCount, sizeof(pair), compareParSecondDec);
	//topZ = (opt->top_Z > distCount) ? distCount : opt->top_Z;

	//FIXME at this version 2014-1-27, topZ always equals to 2
	topZ = distCount;
    /*
	if (posFreqDic[0].second
			>= (opt->firstZ) * (opt->good_enough_percent_support)) {
		//Then there is very high probability that posFreqDic[0].first is the right position
		topZ = 1;
	}
	*/

	if (posFreqDic[0].second < 3 || distCount < 2) //FIXME threshold may change, here we use 1
	{
		lowScoreHandle(otherReadsFiles.postProcessFile, readWhole);
		fprintf(stdout, "\n");
		return 0;
	}

	for (i = 0; i < topZ; i++) {
		uint32_t _shifts[ARRAY_LEN], _shiftCount = 0, j;
		posMatchedRegion[i].pos = posFreqDic[i].first;
		//process posMatchedRegion for that pos
		//if( posFreqDic[i].second > 1)
		{
			uint64_t _pos, k;
			_pos = posFreqDic[i].first;
			for (k = 0; k < dkmerN; k++) { //FIXME Big waist! that process can be combined with previous dictionary process
				if (_pos == positions_shifts[k].position) {
					_shifts[_shiftCount++] = positions_shifts[k].shift;
				}
			}
		}
		if (_shiftCount) //Need sorted,  but shift is actually sorted
		{
			//chaining
			if(old_chaining)
			{
				int start, end, it;
				start = _shifts[0];
				end = start;
				for (it = 1; it < _shiftCount; it++) { //FIXME  still can be optimized. Here we find first longest continues matched region, maybe not the longest
					if (_shifts[it] - end <= KMER_LEN) {
						end = _shifts[it];
					}
				}
				end += KMER_LEN;
				posMatchedRegion[i].region.first = start;
				posMatchedRegion[i].region.second = end;
			}
			else{ //FIXME  New idea, 2014-1-25.
				posMatchedRegion[i].region.first =  _shifts[0];
				posMatchedRegion[i].region.second = _shifts[_shiftCount-1]+ KMER_LEN;
				//fprintf(stdout,"region %d %d\t",posMatchedRegion[i].region.first,posMatchedRegion[i].region.second);

			}
		}
	}

	//FIXME ONLY for short reads, we do this quick check.
	/*
	if (readLength <= SHORT_READ_LENGTH) {
		int mismatchCount, find = 0, minMismatchPosSign;
		uint32_t minMismatchPos, minMismatch = SEMIBIGNUMBER;

		for (i = 0; i < topZ; i++) {
			target = getTarget(pac, l_ref, posMatchedRegion[i].pos, readLength);
			mismatchCount = hammingDistance(target, read, readLength);
			if (mismatchCount <= (int) (readLength * (opt->ert))) {
				//It is a good position
				if (mismatchCount < minMismatch) {

					minMismatchPos =
							(posMatchedRegion[i].pos > BIGNUMBER) ?
									(posMatchedRegion[i].pos - BIGNUMBER) :
									posMatchedRegion[i].pos;
					minMismatchPosSign =
							(posMatchedRegion[i].pos > BIGNUMBER) ? -1 : 1;
					minMismatch = mismatchCount;
					find = 1;
				}
			}
			free(target);
		} //for

		if (find) {
			if(minMismatchPosSign == -1){
				 minMismatchPos = l_ref - minMismatchPos - readLength;
			}
			// get the chrome position
			bns_cnt_ambi(bns, minMismatchPos, readLength, &chrID);
			minMismatchPos = minMismatchPos - bns->anns[chrID].offset + 1;
			fprintf(stdout, "MISM_pos:\t%d\t%llu\t%d\t%s",
					minMismatchPosSign,minMismatchPos,
					opt->a * readLength - minMismatch * opt->b,bns->anns[chrID].name);
			bwaReadsWrite(otherReadsFiles.bwaReadsFile, readWhole);
			fprintf(stdout, "\n");
			return 0;
		}
	} // if read length short enough
	*/

	//Find best alignment
	//FIXME at this version 2014-1-27, topZ always equals to 2

	for (i = 0; i < topZ; i++) {
		int _score = 0, _score_left = 0, _score_right = 0,left_extend_len, right_extend_len, _score_mid =0;
		uint64_t left_start, left_end, left_refLen, right_start, right_end, right_refLen, left_readLen, right_readLen;
		uint8_t *left_read, *right_read,*mid_read;
		uint32_t _len, notUse, mismatchCount,region_len,left_mismatchCount, right_mismatchCount;

		mid_read = &read[posMatchedRegion[i].region.first];
		mid_target = getTarget(pac, l_ref, posMatchedRegion[i].region.first+ posMatchedRegion[i].pos, region_len);
		mismatchCount = hammingDistance(mid_target, mid_read, region_len);
		region_len= posMatchedRegion[i].region.second - posMatchedRegion[i].region.first;
		_score_mid = (region_len - mismatchCount)*(opt->a)-mismatchCount*(opt->b);
		//fprintf(stdout," region_len %d ",region_len);
	    //fprintf(stdout," miscount: %llu  _score_mid %d ",mismatchCount,_score_mid);

		//Prepare reference left and right piece
		left_extend_len =  (posMatchedRegion[i].region.first > 10) ? opt->extend_len : 2;
		left_start = posMatchedRegion[i].pos - left_extend_len;
		left_end = posMatchedRegion[i].pos + posMatchedRegion[i].region.first;
		left_refLen = left_end - left_start;

		right_extend_len = ((readLength - posMatchedRegion[i].region.second) > 10) ? opt->extend_len : 2;
		right_start = posMatchedRegion[i].pos + posMatchedRegion[i].region.second;
		right_end = posMatchedRegion[i].pos + readLength + (right_extend_len);
		right_refLen = right_end - right_start;

		left_target = getTarget(pac, l_ref, left_start, left_refLen);
		right_target = getTarget(pac, l_ref, right_start, right_refLen);

		//Do quick check
		left_mismatchCount  = hammingDistance(&left_target[left_extend_len ], read, posMatchedRegion[i].region.first);
		right_mismatchCount = hammingDistance(right_target, &read[posMatchedRegion[i].region.second], readLength - posMatchedRegion[i].region.second);

		if ((float)left_mismatchCount/posMatchedRegion[i].region.first < 0.25 && (float)right_mismatchCount/(readLength - posMatchedRegion[i].region.second)<0.25){
			_score_left = opt->a * (posMatchedRegion[i].region.first - left_mismatchCount ) - left_mismatchCount*opt->b;
			_score_right = opt->a *(readLength - posMatchedRegion[i].region.second - right_mismatchCount ) - right_mismatchCount *opt->b;
		}
		//fprintf(stdout," left_%d %llu\t",posMatchedRegion[i].region.first,left_mismatchCount );
		//fprintf(stdout,"right_%d %llu\t",readLength - posMatchedRegion[i].region.second,right_mismatchCount );
		else{
			left_readLen =  posMatchedRegion[i].region.first + KMER_LEN/3; //FIXME potential problem here
			right_readLen = readLength - posMatchedRegion[i].region.second + KMER_LEN / 3;

			left_read = &read[0];
			right_read = &(read[posMatchedRegion[i].region.second - KMER_LEN /3]);

			_score_left = dmap_extend(&_len, left_target, left_read, &(parameter.AP), 1, left_refLen, left_readLen, opt);
			_score_right = dmap_extend(&notUse, right_target, right_read, &(parameter.AP), 1, right_refLen, right_readLen, opt);
			//_score = _score_left + _score_right + (posMatchedRegion[i].region.second - posMatchedRegion[i].region.first - (KMER_LEN / 3) * 2) * (opt->a);
		}
		_score = _score_left + _score_right + _score_mid;
		scores[i] = _score;
		printPos(posMatchedRegion[i].pos,_score, readLength , l_ref,bns,1,posFreqDic[i].second);



		//Now we get score for this alignment
		if ( _score > maxScore) {
			maxScore = _score;
		}

		free(left_target);
		free(right_target);
	} //for topZ pos

	//FIXME at this version 2014-1-27, topZ always equals to 2
	if( scores[0] == maxScore ){
		if (posFreqDic[0].second > 4) c1 = 1;
		//fprintf(stdout, " posFreqDic %llu  %d %d ",posFreqDic[0].second,scores[0],scores[1]);
		mapQuality = (float)(c1*250*(scores[0]- scores[1])/scores[0]);
		//printPos(posMatchedRegion[0].pos,maxScore,readLength , l_ref,bns,mapQuality,posFreqDic[0].second);
	}
	else{
		if (posFreqDic[1].second > 4) c1 = 1;
		mapQuality = (float)(c1*250*(scores[1]- scores[0])/scores[1]);
		//fprintf(stdout, " posFreqDic %llu  %d %d ",posFreqDic[1].second,scores[1],scores[0]);
		//printPos(posMatchedRegion[1].pos,maxScore,readLength , l_ref,bns,mapQuality,posFreqDic[1].second);

	}

	fprintf(stdout, "\n");
	return 0;
}

int dmap_map(const char *refFile, const char *readsFile, opt_t *op) {
	int ret, i;
	table_t *kmer_position_table;
	char *tableFn;
	uint8_t *kmerStrandIndicator;
	char *indicatorArrayFn;
	uint8_t *pac;
	bntseq_t *bns;

	uint32_t total = 0, n_seqs = 0;
	uint32_t readsBp = 0;
	seqio_t *seq;
	read_t *readsChunk, *p;
	FILE *bwaReadsFile, *noDkmerReadsFile, *lowScoreReadsFile;
	misReadsFiles otherReadsFiles;
	char *prefix = refFile;
	char *hardMapReadsFn, *noDkmerReadsFn, *lowScoreReadsFn,*postProcessReadsFn;
	clock_t t, mapStart;

	fprintf(stderr, "[dmap_aln] start \n");
	t = clock();
	fprintf(stderr, "[dmap_aln] retrieving index info ...\n");

	//hardMapReadsFn = (char*) calloc(strlen(readsFile) + 30, 1);
	//noDkmerReadsFn = (char*) calloc(strlen(readsFile) + 30, 1);
	//lowScoreReadsFn = (char*) calloc(strlen(readsFile) + 30, 1);
	//strcpy(hardMapReadsFn, "easyMapReads_");
	//strcat(hardMapReadsFn, readsFile);
	//strcpy(noDkmerReadsFn, "noDkmerMapReads_");
	//strcat(noDkmerReadsFn, readsFile);
	//strcpy(lowScoreReadsFn, "lowScoreMapReads_");
	//strcat(lowScoreReadsFn, readsFile);

	postProcessReadsFn =(char*) calloc(strlen(readsFile) + 30, 1);
	strcpy(postProcessReadsFn, "postProcessReadsFn_");
	strcat(postProcessReadsFn, readsFile);

	//otherReadsFiles.bwaReadsFile = fopen(hardMapReadsFn, "w");
	////otherReadsFiles.noDkmerReadsFile = fopen(noDkmerReadsFn, "w");
	//otherReadsFiles.lowScoreReadsFile = fopen(lowScoreReadsFn, "w");
	otherReadsFiles.postProcessFile = fopen(postProcessReadsFn, "w");

	tableFn = (char*) calloc(strlen(refFile) + 10, 1);
	strcpy(tableFn, refFile);
	strcat(tableFn, ".tb");
	indicatorArrayFn = (char*) calloc(strlen(refFile) + 10, 1);
	strcpy(indicatorArrayFn, refFile);
	strcat(indicatorArrayFn, ".idc");

	bns = bns_restore(prefix);
	//When we do mapping, we only need + strand reference

	pac = calloc((bns->l_pac / 2) / 4 + 1, 1);

	if (pac == 0) {
		fprintf(stderr, "[dmap_aln] insufficient memory!\n");
		exit(1);
	}
	fread(pac, 1, (bns->l_pac / 2) / 4 + 1, bns->fp_pac);
	fprintf(stderr, "[dmap_aln] number of characters of reference  %llu \n",
			bns->l_pac);

	kmer_position_table = table_restore(tableFn);
	kmerStrandIndicator = indicatorArray_restore(indicatorArrayFn);
	fprintf(stderr, "[dmap_aln] retrieve success! %.2f sec\n",
			(float) (clock() - t) / CLOCKS_PER_SEC);

	mapStart = clock();


	seq = seq_open(readsFile);
	/*Core loop, reads mapping. To save memory, sequentially process chunk of
	 * reads then write results to disk
	 * */
	readsChunk = (read_t*) calloc((op->chunck_size) + 1, sizeof(read_t));

	while (kseq_read(seq->ks) >= 0) {
		total += 1;
		p = &readsChunk[n_seqs++];

		{ //Read a sequence
			p->name = strdup((const char*) seq->ks->name.s);
			p->len = seq->ks->seq.l;
			readsBp += p->len;

			//copy read sequence
			p->seq = (ubyte_t*) calloc(p->len, 1);
			for (i = 0; i != p->len; i++)
				p->seq[i] = nst_nt4_table[(int) seq->ks->seq.s[i]];

			// copy read quality
			if (seq->ks->qual.l) {
				p->qual = (ubyte_t*) strdup((char*) seq->ks->qual.s);
				p->qual_len = seq->ks->qual.l;
			}
		}
		//fprintf(stdout, "\n");
		//fflush(stdout);
		if (n_seqs == (op->chunck_size) || readsBp > (op->chunck_size) * 100) {

			fprintf(stderr, "[dmap_aln] read %d sequences ...", total);
			t = clock();
			ret = dmap_map_core(bns, kmer_position_table, kmerStrandIndicator,
					(bns->l_pac) / 2, pac, readsChunk, n_seqs, op,
					otherReadsFiles);
			fprintf(stderr, "%.2f sec\n",
					(float) (clock() - t) / CLOCKS_PER_SEC);
			n_seqs = 0;
			readsBp = 0;
			readsChunk = (read_t*) calloc((op->chunck_size) + 1,
					sizeof(read_t));
			fprintf(stderr, "unableMapReadsCount: %llu \n",
					unableMapReadsCount);
			fprintf(stderr, "scoreLowThresCount: %llu \n", scoreLowThresCount);
		}
	}
	fprintf(stderr, "[dmap_aln] read %d sequences ...\n", total);
	ret = dmap_map_core(bns, kmer_position_table, kmerStrandIndicator,
			(bns->l_pac) / 2, pac, readsChunk, n_seqs, op, otherReadsFiles);

	//DEBUG ONLY
	fprintf(stderr, "\nNinReadsCount: %llu \n", NinReadsCount);
	fprintf(stderr, "TOTAL READS: %llu \n", total);
	//fprintf(stderr, "unableMapReadsCount: %llu \n", unableMapReadsCount);
	fprintf(stderr, "scoreLowThresCount: %llu \n", scoreLowThresCount);
	fprintf(stderr, "notHasDKmer: %llu \n", notHasDKmer);

	fprintf(stderr, "[dmap_aln] map time %.2f sec\n",
				(float) (clock() - mapStart) / CLOCKS_PER_SEC);
	//DEBUG ONLY
	//fclose(otherReadsFiles.bwaReadsFile);
	//fclose(otherReadsFiles.lowScoreReadsFile);
	//fclose(otherReadsFiles.noDkmerReadsFile);
	fclose(otherReadsFiles.postProcessFile);

	free(tableFn);
	free(indicatorArrayFn);
	free(pac);
	free(kmer_position_table);
	free(kmerStrandIndicator);
	free(bns);
	free(seq);
	return 0;
}

static void generateNeighbors(const uint32_t number,uint32_t *neighbors,uint32_t maskALL)
{   int i;
	for(i = 0; i < 4; i++){
		neighbors[i] = number + i;
		neighbors[i] &= maskALL;
	}
}
int dmap_map_core(bntseq_t *bns, table_t *kmer_position_table,
		uint8_t *kmerStrandIndicator, const int64_t l_ref, const uint8_t *pac,
		read_t *readsChunk, uint32_t n_seqs, opt_t *op,
		misReadsFiles otherReadsFiles) {
	uint32_t i = 0;
	int k, thatReadHasDistKmer,loop_p;
	read_t *p, *temp;
	uint32_t firstKey = 0, oldKey, nowKey;
	uint32_t maskALL = (1 << (KMER_LEN * 2)) - 1;


	uint64_t shifts[ARRAY_LEN]; //Positions that has distinct Kmer
	uint64_t positions[ARRAY_LEN]; //Kmer position in reference
	pos_shift positions_shifts[ARRAY_LEN];

	uint32_t dkmerN = 0; //number of distinct kmer.  Always < =  opt->firstZ
	position_t position;
	twoPositions_t tempTwoPos;
	int jumpCount = 0;

	uint32_t neighbors[4], checkKmer;

	//initialize parameter
	alnParam parameter;
	parameter.allowMismatch = op->allowMismatch;

	if (KMER_LEN == 16)
		maskALL = 0xffffffff;


	while ((uint32_t) i < (uint32_t) n_seqs) {

		thatReadHasDistKmer = 0;
		p = &readsChunk[i++]; //A read.

		//Find first distinct kmer in that read
		//if ( i % 1000 == 0 ) fprintf(stderr,"reads: %llu \n",i);
		for (k = 0; k < KMER_LEN; k++) {
			firstKey <<= 2;
			firstKey += (uint32_t) (tripN(p->seq[k])); //FIXME, N is replaced randomly by A C T G
		}
		oldKey = firstKey;

		//FIXME The first kmer not handled properly
		//FIXME  firstZ need be handled more properly
		//op->firstZ = (p->len) * (op->distinctPercent + op->ert);
		//op->firstZ  = 20;
		//Get all the distinct kmer in that read
		for (k = KMER_LEN; k < p->len; k++) {

			oldKey <<= 2;
			nowKey = (oldKey + (uint32_t) tripN(p->seq[k]));
			nowKey &= maskALL;
			generateNeighbors(oldKey,&neighbors, maskALL);
			oldKey = nowKey;
			//fprintf(stderr, " now: %s  \n",int2bin(nowKey, NULL));

			//for(loop_p = 0; loop_p < 4; loop_p++)
			//	fprintf(stderr, " nei: %s  \n",int2bin(neighbors[loop_p], NULL));

			//fprintf(stderr, "\n");
			//FIXME idea change 2014-Feb-4, add a for loop, len 4
			//for(loop_p = 0; loop_p < 4; loop_p++){
				//checkKmer = neighbors[loop_p];
				checkKmer = nowKey;
				//FIXME I need to introduce jump here . 2013-11-22
				if (isUniqueArray1(kmer_position_table, checkKmer)) {

					shifts[dkmerN] = (k - KMER_LEN + 1);
					position = getValueArray1(kmer_position_table, checkKmer);

					//positions[dkmerN] = position - shifts[dkmerN];
					positions_shifts[dkmerN].position = position - shifts[dkmerN];
					positions_shifts[dkmerN].shift = shifts[dkmerN];

					if (getKmerStrandIndicatorBit(checkKmer, kmerStrandIndicator) > 0) //that distinct keymer located in - strand
							{
						//positions[dkmerN] = position - shifts[dkmerN] + BIGNUMBER;
						positions_shifts[dkmerN].position = position
								- shifts[dkmerN] + BIGNUMBER;
						positions_shifts[dkmerN].shift = shifts[dkmerN];

					} else {
						//positions[dkmerN] = position - shifts[dkmerN];
						positions_shifts[dkmerN].position = position
								- shifts[dkmerN];
						positions_shifts[dkmerN].shift = shifts[dkmerN];

					}
					dkmerN++;
					thatReadHasDistKmer = 1;
				} //if one kmer is distinct
			//}

			//if (dkmerN > op->firstZ) { //FIXME may need break earlier
			//	break; //break for
			//}

		} //for, all kmer in that read

		if (thatReadHasDistKmer == 0) {
			//pass it to unmapped
			notHasDKmer++;
			add2UnmapedReads(otherReadsFiles.postProcessFile, p);
		} else {
			fprintf(stdout, "%s ", p->name);
			Aln(bns, parameter, &positions_shifts, dkmerN, l_ref, pac,
					p->len, p->seq, op, otherReadsFiles, p);
		}
		dkmerN = 0;
		//exit(1);
	} //while
	free_read_seq(n_seqs, readsChunk);
	return 0;
}

